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ABSTRACT 

LI  phase  measurements  collected  by  a  miniaturized  and  low  cost  GPS  logging  device  are  the  basis  for  a 
navigation  approach  aiming  to  precisely  measure  (flight)  trajectories  independently  from  any  auxiliary 
system.  Forming  single  differences  between  two  measurements  taken  by  one  moving  receiver  at  two 
different  times  allows  achieving  relative  precision  in  the  low  decimeter  range  over  time  intervals  of  up  to 
several  minutes.  Neither  space  nor  ground  based  augmentation  systems,  a  second,  nearby  base  receiver  or 
any  (static)  initialization  patterns  are  required  by  the  method.  This  fact  significantly  reduces  the 
complexity  to  be  handled  by  the  user  often  operating  under  adverse  field  conditions  even  at  extremely 
remote  locations.  The  approach  is  taking  advantage  of  canceling  ambiguities  instead  of  striving  for 
estimating  these  unknown  quantities  even >  phase  measurement  is  biased  with.  The  constitutive  navigation 
equations  are  derived  and  theoretic  aspects  of  various  error  sources  restricting  possible  processing 
inten’als  are  discussed  in  this  paper.  The  method  is  validated  by  both  static  and  dynamic  reference  data. 
Finally  first  results  from  a  GPS  campaign  at  Kerguelen  Archipelago  to  measure  the  dynamic  soaring  of 
Albatrosses  using  the  time-difference  method  are  presented  and  exemplary  data  are  processed  using  the 
time-difference  method. 


1.0  INTRODUCTION 

The  context  the  present  navigation  approach  is  developed  in  is  a  rather  unconventional  task:  Exploring  the 
question  how  (Wandering)  Albatrosses,  the  world’s  largest  birds,  can  fly  thousands  of  miles  over  open 
seas  without  flapping  their  wings  [21].  Considering  the  absence  of  vertical  thermal  up  winds  as  utilized  by 
birds  (or  glider  planes)  soaring  over  solid  land  mass,  this  capability  becomes  even  more  fascinating. 
Extensive  research  has  been  conducted  on  the  theoretical  aspects  of  the  birds’  flight  technique  named 
“dynamic  soaring”  -  a  cyclic  maneuver  performed  in  the  shear  wind  field  right  above  the  water  surface 
[13].  The  maneuver  and  the  geometry  of  an  individual  cycle  are  illustrated  in  Figure  1.  Measuring  suchlike 
maneuvers  is  complicated  by  multiple  circumstances.  For  flight  mechanical  analysis  position  fixes  are 
required  with  a  relative  precision  to  the  maneuver’s  starting  point  in  the  decimeter  range.  High  dynamics 
require  high  sampling  rates  of  at  least  10  Hz.  The  sensor  used  must  operate  completely  autonomously  and 
has  to  be  lightweight  (<  100  g)  and  small  (<  150  x  50  x  20  mm3).  The  weight  restriction  entails  low  power 
consumption  for  a  runtime  of  up  to  50  -  70  hours  just  as  budget  limitations  and  the  risk  of  not  recovering 
the  device  imply  the  use  of  low  cost  sensors.  The  latter  requirements  are  readily  met  by  standard  single 
frequency  GPS  receiver  modules  as  used  in  car  navigation  systems  or  latest  generation  cell  phones. 
However,  sufficient  precision  cannot  be  provided  and  satellite  or  ground  based  augmentation  systems 
cannot  be  addressed  due  to  the  extremely  remote  habitat  of  the  Albatross,  compare  Figure  2.  Some  of  the 
GPS  LI  receiver  modules  output  raw  data,  i.e.  pseudorange,  carrier  phase,  Doppler  frequency  and  raw 
signal  strength  or  SNR  values,  with  up  to  10  Hz.  Logging  these  raw  (phase)  measurements  for  sufficiently 
precise  postprocessing  instead  of  sticking  to  the  code  based  online  navigation  solution  qualifies  these 
modules  as  appropriate  sensors  for  the  described  flight  measurement  task.  Note  that  the  online  navigation 
solution  is  restricted  to  sampling  rates  well  below  the  rate  the  raw  measurements  are  provided  with. 
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ORGANIZATION 


Conventional  phase  based  processing  techniques  rely  on  the  presence  of  data  collected  by  a  nearby  (<  10 
km)  base  station  and  (mostly  static)  initialization  patterns  of  the  roving  receiver.  Both  cannot  be  provided 
in  the  scope  of  the  current  applications:  birds  neither  remain  static  nor  close  to  a  base  station  and  can 
perform  maneuvers  inevitably  causing  signal  shadowing  resulting  in  complete  loss  of  lock  undoing  any 
kind  of  initialization.  Hence  an  alternative  way  to  process  the  LI  phase  measurements  with  no  need  for  a 
base  station  or  initialization  patterns  is  discussed  in  this  paper.  The  price  to  be  paid  for  these  advantages  is 
a  limitation  in  processing  time  to  a  view  minutes.  A  detailed  description  of  the  theoretic  concept  is  given 
as  well  as  a  practical  validation  by  static  experiments  and  flight  tests  with  an  experimental  aircraft.  Finally 
results  from  our  Winter  2008/09  Albatross  measurement  campaign  to  Kerguelen  Archipelago  are 
presented. 


(a)  Flight  over  water  surface 


(b)  Geometry  of  a  cycle,  flown  within  approx.  7 
- 10  sec 


Figure  1:  Dynamic  soaring  of  Albatrosses  [13] 


Figure  2:  The  habitat  of  the  Wandering  Albatross,  Kerguelen  Islands  (red  circle), 

with  SBAS  availabilitiy  [6] 
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2.0  THEORETICAL  BACKGROUND 

2.1  The  Observable:  LI  Phase  Ranges 

The  present  method  is  based  on  the  L 1  phase  range  observable  generated  in  the  phase  lock  loop  (PLL)  of  a 
GPS  receiver.  The  nominal  (constant)  frequency  /  of  the  LI  carrier  wave  is  1575.42  MHz.  Let  cps (t) 

[cycles]  designate  the  phase  of  the  satellite  emitted  wave  as  perceived  by  the  receiver  and  (pR(t )  [cycles] 
denote  the  phase  of  the  (not  yet  Doppler  compensated)  receiver  replica  signal.  With  t  [s]  as  an  epoch  in 
GPS  system  time  reckoned  from  an  initial  epoch  t0  (=  0)  when  acquiring  lock  to  the  respective  PRN  one 
can  state: 


fP^  Po(0 

c 

<ps0{t)  =  -Ss{t)f+<ps{t0 ) 

(1) 

<PR(t)  =  ft-<Po(t) 

cpR(t)=-S"(t)f  +  (pK(tQ) 

(2) 

Here,  tpg  and  <pR  are  the  satellite  and  receiver  initial  phase  biases  at  tQ  afflicted  with  the  respective  clock 

errors  Ss  and  SR  [s]  ( 5  >  0  «->  clock  reading  ahead  of  GPS  system  time).  The  symbol  p  [m]  is  the 
geometric  distance  between  satellite  and  receiver;  c  [m/s]  denotes  the  nominal  signal  propagation  speed 
in  vacuum.  The  beat  phase  </>RS  [cycles]  is  now  given  by 

c 

Note  that  due  to  the  initial  distance  between  satellite  and  receiver  p(t0),  Eq.  (3)  theoretically  yields  a  large 

decimal  number  for  the  beat  phase  </>RS  (t).  However  when  acquiring  lock  to  the  respective  PRN  only  the 
fractional  part  of  this  number  can  be  measured,  the  initial  integer  number  N  of  cycles  between  satellite 
and  receiver  is  unknown.  Note  that  this  ambiguity  is  not  time  dependent  as  long  as  the  phase  is  locked 
continuously:  N^N(t).  Hence,  to  model  the  numerical  value  actually  output  by  the  receiver,  this 
unknown  but  constant  integer  term  N  has  to  be  subtracted  from  the  right-hand  side  of  Eq.  (3).  Moreover,  to 
obtain  a  value  directly  corresponding  to  a  geometric  range,  the  result  is  multiplied  by  -2,  (the  LI 
wavelength)  yielding  the  (pseudo)  phaserange  ®  [m]: 


^)(Q-/7(0  +  c(^j?(0~^5(0)+ a)  (4) 

~N] 

Neither  the  initial  satellite  nor  receiver  phase  <ps(t0),  (pR{t0)  are  known.  Hence  the  non- integer  ambiguity 
term  N'  is  introduced  as  shown  in  Eq.  (4)  to  finally  rewrite  for  the  LI  phaserange: 

®(f)  =  p(t)+c(dR{t)-5s(tj)+\N'  (5) 

In  literature  the  non-integer  ambiguity  term  N'  is  frequently  replaced  by  the  integer  ambiguity  N .  This  is 
not  precise  in  a  strict  sense  but  the  fractional  receiver  and  satellite  phase  terms  only  drop  out  and  a  truly 
integer  ambiguity  term  is  obtained  when  forming  double  differences.  However,  the  non-integer  character 
of  the  non-differenced  phase  observations  does  play  an  important  role  for  procedures  such  as  precise  point 
positioning  [7].  Keep  in  mind  that  the  phase  observable  changes  in  the  same  sense  as  the  C/A  code 
pseudorange  (negative  Doppler)  but  is  significantly  less  noisy. 
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2.2  Basic  Concept 

The  principal  idea  of  the  present  approach  is  the  fact  that  the  ambiguity  N'  is  time-invariant,  provided 
lock  of  phase.  Hence,  differencing  phase  observations  to  the  same  satellite  S  across  two  epochs  tb  and  tt 
allows  for  canceling  N1 : 


biD O  =  O,.  -  0A  =  biDp  +  cbiD8R  +  /L,  biDN'  (6) 

o 

with  D  denoting  single  temporal  differences.  The  right  hand  side  subscript  i  shall  be  a  short  equivalent 
for  the  time  argument  (/;.)  during  the  following  derivations.  Note  that  the  satellite  clock  error  DSs  is 
neglected  in  Eq.  (6). 

The  shown  way  to  cancel  ambiguities  is  the  background  for  a  diversity  of  applications.  In  tightly  coupled 
GPS/INS  systems  time  differenced  double  differences  (so  called  triple  differences  across  two  receivers, 
two  satellites  and  two  epochs)  can  support  the  dynamics  estimation  for  attitude  computation  [5].  In  a 
similar  context,  also  carrier  phases  directly  differenced  between  subsequent  epochs  can  be  used  instead  of 
the  noisier  delta-range  (Doppler)  measurements  to  improve  velocity  and  attitude  information  without  the 
need  for  a  base  station  [20],  Triple  differences  can  also  be  used  for  carrier  phase  cycle  slip  detection  [10]. 
Further,  precise  baseline  computation  between  a  base  and  a  roving  receiver  provided  that  there  are  at  least 
7  satellites  (if  only  phase  data  are  used)  in  view  is  possible  with  triple  differences  [8].  Time-differences 
used  in  a  stand-alone  GPS  application  are  used  to  process  static  data  for  gun-laying  applications  [17].  This 
approach,  enhanced  by  a  loop  misclosure  procedure,  can  be  applied  to  static  measurements  from  civil 
receivers  since  selective  availability  has  been  switched  off  [1]. 

Using  time  differences  for  processing  kinematic  data  is  an  unconventional  approach  emerged  from  the 
need  for  a  high  quality  but  low  effort  navigation  solution.  The  navigation  equations  are  derived  in  the 
following. 


2.3  Navigation  Equations 

Assume  the  position  x  of  the  receiver  to  be  known  at  a  (base)  epoch  tb  .  Then  one  can  write  for  the  base 
vector  pointing  from  the  position  at  the  base  epoch  to  the  position  at  the  epoch  of  interest  tt 


bbi 


(7) 


Figure  3  illustrates  the  geometry  of  the  problem. 

Besides  the  position,  the  receiver  clock  error  with  respect  to  GPS  system  time  has  to  be  determined, 
yielding  the  combined  position  and  time  (PT)  variable: 

^  =  [xt,cSr]T  (8) 
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Here  the  clock  error  is  scaled  to  range  by  multiplication  with  the  signal  propagation  speed  for 
convenience.  Now  one  can  set  up  the  enhanced  base  vector: 

r=^-^b  (9) 


In  order  to  determine  |1' ,  time  differences  hl  Dd)  between  two  phase  observations  at  tb  and  /,  as 
introduced  in  Eq.  (6)  are  used.  For  precise  navigation,  the  phaserange,  Eq.  (5),  has  to  be  enhanced  by 
models  for  signal  delays  caused  by  ionospheric  and  tropospheric  refraction,  /  and  T : 


®,  =  pt  +  cdf  +  ^ N'+Tf  -  /, 


(10) 


Here  the  symbol  A  is  introduced  to  distinguish  modeled  from  measured  values  denoted  ~  in  the  remainder. 
The  geometric  range  is  calculated  by 


A  = 


(11) 


where  X,  is  the  satellite  position  at  the  time  of  signal  emission  corresponding  to  the  epoch  li .  This 

position  (in  conjunction  with  the  respective  emission  time)  is  extracted  from  the  Ephemeris  data  of  the 
respective  satellite.  Now  the  time  differenced  phase  observations  can  be  rewritten  as 


/XI)  =  blDp  +  c  " D8r  -hi  DI  +bl  DT  =  f\£t , X, ,  t, 


>  ’  X b  5  tb  ) 


(12) 


with  as  the  only  unknown.  Note  that  in  a  strict  sense  D<b  is  not  a  function  of  p  but  of  both  the  PT 
solution  at  the  current  epoch  and  the  (known)  base  epoch.  In  order  to  solve  for  at  least  four 
measurements  are  required: 


i'DO=wDO(^.) 


(13) 


With 


bi D<*>  =  [bi Dab  bi  m2  ...  biD<t>m]T ;m>4 


RTO-EN-SET-1 16(2010) 


4-5 


Conceptual  Approach  for  Precise  Relative  Positioning 
with  Miniaturized  GPS  Loggers  and  Experimental  Results 


ORGANIZATION 


Just  as  for  standard  single  point  processing,  this  (over-determined)  set  of  equations  is  solved  by  the 
nonlinear  least  squares  method.  For  that  purpose,  the  right  hand  side  of  Eq.  (13)  has  to  be  linearized: 


-■z>6fo,0=  W0®fe,o,O+  H 


f/,0 


AS, 


(14) 


The  linearization  point  0  is  either  the  result  of  the  last  iteration  cycle  of  the  current  epoch  or  the  final 
result  of  the  previous  epoch.  The  Jacobian  writes  elaborately  as 


H 


dbiD<S> 


\i,  0 


f,.o 


With  Eqs.  (10)  and  (12)  the  following  relationships  hold 

5h:D(b  ~  dp,  dbiD<& 


fiX,.  fiX; 


:  -e. 


4^) 


-i 


S^DO1 

8h,D<b] 

dbiD& 

8lnD(b' 

fix,- 

dhiD<b 2 

5y, 

db,D& 

dz, 

dbiD& 

d(cS«) 

dbiD<i>2 

fix,- 

dy, 

dzj 

died}) 

dh,D<bm 

dbiD&n 

dbiDi hm 

8b‘D(bm 

fix,- 

fiy, 

fiz,. 

died}) 

(15) 


(16) 


where  e,  is  the  unit  vector  pointing  from  the  rover  position  at  t,  to  the  satellite,  compare  Figure  3.  Spatial 
variations  of  the  atmospheric  models  are  neglected  when  linearizing  as  these  effects  are  compensated  for 
by  the  iteration.  With  Eq.  (16)  the  Jacobian  rewrites  to 


H 


\i.  0 


(17) 


Now  one  obtains  with  least  squares  for 

and  iteration 

i,k+ 1  —  ^i,k  1  /~xs3i 

finally  yields  %, .  Note  that  even  so  this  solution  directly  drops  out  of  the  solution  process,  is  precise 
relative  to  %h  only  and  its  absolute  accuracy  depends  on  the  one  of  \b  itself.  Consequently,  \\h‘  is  the 
intrinsic  solution  of  the  problem.  This  will  become  even  more  evident  in  Eq.  (28). 

The  shown  solving  procedure  minimizes  the  residuals  fj  between  measured  and  modeled  observations  in 
a  least  squares  sense: 


(18) 


.  A? 


(19) 
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£(/•'?=  min  <2°) 

7=1 

f  i  = bi  J  - bi  £><b j  1 

If  the  inevitable  experimental  measurement  errors  are  uncorrelated,  have  a  mean  of  zero  and  a  constant 
variance,  the  least  squares  estimate  of  has  the  minimum  variance  of  all  estimates  that  are  linear 
combinations  of  the  observations.  In  this  sense  it  is  the  optimal  solution. 


2.4  Error  Analysis 

As  indicated  in  Eq.  (10)  the  phase  range  model  is  enhanced  by  terms  taking  into  account  atmospheric 
signal  propagation  delays.  However,  these  effects  cannot  be  eliminated  completely.  The  remaining  effects 
s  compose,  together  with  other  contributions,  the  non-modeled  range  error  % : 


cly  =6i+j,  (22) 

Xi  =  —cSf  +  Ei  +  st  +  noise  +  multipath 

This  error  is  made  up  by  the  scaled  satellite  clock  bias  5s ,  the  satellite  position  offset  due  to  uncertainties 
in  the  satellite  ephemeris  E ,  the  described  remaining  atmospheric  delays  s  ,  measurement  noise  and 
multipath.  When  forming  time  differences  this  yields 

hlD$>  =  b,D®  +  b,Dx  (23) 

biDX  =  -c  hiD5s  +  hiDE  +  biDs  +  biD  noise  +bi  D  multipath 

Except  for  the  noise  ( <rVnoise  «  <rnoise )  and  multipath,  the  range  error  components  are  highly  temporally 
correlated  which  allows  the  following  “linearization”: 


b‘DZ  = 


dSs 

dE 

de] 

/ 

\  ds 

- c - 

+  — 

+  — 

it: 

t, l  )  +  p 

dt 

dt 

dt 

h 

O  /  ^  A 

O  X 

b;  +h'D  noise  +hlD  multipath 

h 


(24) 


The  atmospheric  error  stemming  from  spatial  variations  occurs  in  any  kind  of  differential  processing  and 
remains  small  with  short  baseline  length.  As  a  matter  of  fact,  the  raise  of  the  noise  level  cannot  be  avoided 

but  is  smaller  by  a  factor  of  1/V2  than  when  working  with  double  differences.  However,  the  drift  terms 
are  the  bottleneck  problem  for  the  time  difference  approach.  In  order  to  keep  this  effect  small,  precise 
correction  data  are  utilized.  These  data  not  only  decrease  the  absolute  range  error  but  also  increase  its 
temporal  correlation.  Usually  only  little  attention  is  paid  to  this  fact  which  is  key  for  extending  the  precise 
processing  time  when  working  with  time  differences  [15].  However,  there  is  yet  another  error  source 
affecting  the  relative  precision  of  the  solution  -  the  bias  in  from  the  true  position  and  time.  In  order  to 

analyze  this  effect,  Eq.  (14)  is  linearized  in  both  0  and  the  true  base  position  and  timely  0 : 


iZj6feI,§J  =  wZ)6feI,0>^>0)+H5(oA§I.-H54 


(25) 


with  the  Jacobian  according  to  Eq.  (17)  and  A^h  =  %b-£sb0.  Referring  to  Eqs.  (10),  (12)  and  (17) 
one  gets  for  the  second  Jacobian 
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= 


8biV  6fe(.,§4) 


3$b 


-e 


b,  0 


T 

-em  1 


(26) 


The  minus  sign  is  used  for  convenience  only  in  order  to  obtain  analogy  between  H-  and  H  ^  n .  Figure  3 
illustrates  the  geometry  of  the  unit  vectors.  The  matrices  1 1 ,  and  FF  are  closely  related  which  allows 
again  to  linearize: 


HV.=Hi,+"d'» -'.)=HVl  + 


-e 


i,0 


-e 


i,0 


(h-h) 


(27) 


Note  that  Eq.  (27)  only  takes  into  account  changes  due  to  the  time  elapsed  between  base  and  rover  epoch 
not  due  to  the  motion  of  the  receiver.  The  former  does  have  a  significant  impact  on  the  unit  vectors  due  to 
the  motion  of  the  satellites  themselves  whereas  the  latter  is  neglected  due  to  the  far  distance  between 
receiver  and  satellites.  Inserting  Eq.  (27)  in  Eq.  (25)  yields 

wv6feI.,^)=wv6feI.0,§M)+ 

tIFj.U.  A^).H4(t,  th).\Zh  (28) 

AP6'  b‘DZgeo 

The  last  term  of  Eq.  (28)  is  unknown  and  can,  in  analogy  to  Dx ,  Eq.  (23),  be  interpreted  as  a  geometric 
range  error  caused  by  variable  satellite  geometry: 


bi 


biDZl=-e{TA^b{ti-tb) 

=[biDxl,...,biDXt  Y 


(29) 


Just  as  the  non-modeled  range  error  Dx ,  the  geometric  range  error  causes  an  error  in  the  relative 
navigation  solution. 


Figure  4  aims  to  illustrate  this  effect:  If  the  position  at  tb  is  afflicted  by  an  offset  S\h ,  the  vector  Iv"  is 
first  of  all  translated  by  Sx.h  as  well.  This  translation  causes  a  shift  of  x;  by  S\h  which  is  no  error  as 
relative  precision  is  yet  maintained.  Flowever,  just  as  every  other  range  error,  also  V/_gco  causes  a 

distortion  of  l/'  to  b/>'*,  indicated  by  <5b/>'  in  the  figure.  This  does  degrade  relative  precision.  It  is 
interesting  to  note  that  the  variation  of  the  satellite  geometry  over  time  has  virtually  no  impact  on  the 
relative  solution  if  the  initial  position  and  time  is  known  exactly,  §^b  =  0 ,  no  matter  how  long  time 
intervals  are  processed. 
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Figure  4:  Impact  of  a  bias  in  the  initial  epoch’s  position  on  relative  precision 


In  order  to  estimate  the  impact  of  the  various  range  errors  on  p  1  ,  Eq.  (23)  is  enhanced  to 


D&  =bi  /)([) '  +  °‘DzJ  +m  D%, 


j  _i_  bi  , 


geo 


'*>zL 


(30) 


The  total  range  error  is  a  superposition  of  different  random  processes.  Each  of  them  has  a  different  degree 
of  temporal  autocorrelation.  For  instance  the  change  of  the  error  caused  by  imprecise  ephemeris  is  much 
slower  than  the  change  of  the  measurement  noise.  As  only  short  time  intervals  tt  -  tb  are  processed,  the 
slowly  changing  error  contributions  can  be  considered  as  systematic  but  unknown  biases  in  the 
measurements.  In  other  terms,  the  individual  range  errors  can  be  assumed  to  be  non-ergodic 

random  processes  with  constant  and  equal  variance  and  statistic  distribution  but  a  different,  unknown 
expectation  which  is  growing  with  time  for  each  measurement.  This  fact  significantly  complicates  the 

estimation  of  the  error  in  P  .  In  the  scope  of  this  paper,  the  quality  of  the  solution  is  validated  via 
comparison  with  reference  solutions.  In  addition,  a  simplistic  estimation  of  the  error  propagation  is 
provided  by  referring  to  the  concept  of  Dilution  of  Precision  (DoP)  as  known  from  standard  single  point 
processing:  Ignoring  all  systematic  errors  in  the  measurements  for  the  lack  of  better  knowledge,  a  rough 
estimation  of  the  covariance  matrix  of  the  least  squares  solution  is  obtained  by 


(31) 


For  the  number  of  measurements  exceeding  the  number  of  unknowns,  m»n  =  4,  the  variance  of  the 
measurements  is  estimated  from  the  residuals  by 


z;./ 

177-4 


(32) 


Using  the  residuals  from  the  solution  between  subsequent  epochs,  this  represents  the  stochastic 

error  component  of  the  measurements.  Using  the  residuals  of  the  over-all  solution,  tb ,  ti ,  one  obtains  a 
quantity  representing  the  error  drift. 


The  DOP  value  referring  to  position  accuracy  is  defined  by 

PDOP  = 


(33) 
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By  comparing  coefficients  an  estimate  for  either  the  3D  noise  component  or  the  3D  error  drift  is  obtained 
by 

=^*'pD°p  <34) 

Note  that  the  Jacobian  when  working  with  time  differences  as  indicated  in  Eq.  (17)  coincides  with  the  one 
used  when  processing  single  points.  Hence  the  numerical  DOP  values  characterizing  the  quality  of  a  single 
point  solution  also  hold  for  the  time  difference  solution. 


3.0  IMPLEMENTATION 

3.1  Two  different  implementation  strategies 

The  above  equations  describe  the  core  algorithm  of  the  navigation  task  -  the  relative  solution  between  two 
individual  epochs.  In  order  to  implement  a  tool  suitable  for  kinematic  trajectory  measurement,  two 
different  strategies  can  now  be  applied  as  shown  in  Figure  5. 


s'  V 

/  ' 

X  S 

/  \ 


SPP 


o 


(a)  Over-all  solution 


(b)  Accumulated  solution 


Figure  5:  Different  implementation  strategies  of  the  time-difference  approach 

In  either  case,  the  initial  solution  qh  is  estimated  by  code-based  single  point  processing  (SPP)  and 
therefore  typically  afflicted  by  an  absolute  error  in  the  meter  range  ( (Kh ).  The  first  strategy  to  obtain  the 

complete  trajectory  is  to  calculate  the  “over-all  solution”  p  1  for  each  epoch  as  illustrated  in  Figure  5(a). 
Starting  from  the  core-algorithm,  this  approach  suggests  itself  -  each  epoch  is  calculated  independently 
from  all  other  observations  and  it  is  straightforward  to  judge  the  quality  of  the  solution  (error  drift)  via  Eq. 

(34)  based  on  the  (raising)  residual  level  ^  -b'Dd>JJ  ,  Eqs.  (21)  and  (32).  However,  there  are 

several  caveats:  The  same  subset  of  satellites  must  be  available  at  both  tb  and  /, .  Considering  that 
satellites  might  be  lost  on  the  way  due  to  shadowing  or  antenna  tilting,  this  is  a  severe  restriction. 
Moreover  signals  which  are  once  afflicted  by  cycle  slips  (r.t.  the  next  section)  must  not  be  reused,  unless 
there  is  a  reliable  repair  algorithm.  Another  problem  can  occur  with  the  drop-out  of  individual  satellites 
which  can  induce  small  but  sudden  jumps  in  the  trajectory  -  what  may  impede  velocity  determination  via 
differentiation  in  the  position  domain.  The  alternative  is  to  accumulate  the  position  increments  obtained  by 

solving  between  subsequent  epochs  P"-1’"  as  demonstrated  by  Figure  5(b).  Accumulating  error  afflicted 
quantities  always  implies  the  risk  of  getting  a  random-walk  like,  unbound  error.  As  demonstrated  by 
various  experiments  and  shown  theoretically  in  [16],  this  is  not  the  case  here  as  the  errors  from  the 
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incremental  solution  perfectly  cancel  when  accumulating  to  the  final  solution.  It  turns  out  that  the  over-all 
and  the  accumulated  solution  are  identical  in  a  first  order  approximation  given  equivalent  constellation 
geometry.  Now  the  various  advantages  of  the  accumulation  strategy  become  evident:  The  common  subset 
of  satellites  between  two  subsequent  epochs  tn_x  and  tn  can  be  used  yielding  generally  improved 
geometries  compared  to  the  subset  available  at  ti  and  tb  (more  common  satellites  available).  Furthermore, 

cycle-slips  do  not  propagate  to  subsequent  epochs  which  makes  a  discrimination  between  simple  outliers 
and  cycle-slips  obsolete  when  monitoring  the  integrity  of  the  solution.  Flowever,  estimating  the  error  drift 
is  based  on  the  residual  level  of  the  over-all  solution  -  a  quantity  which  is  no  longer  directly  available  but 
must  be  reconstructed  by  accumulating  the  residuals  of  the  incremental  solutions.  This  leads  to  overly 
pessimistic  error  drift  estimates.  Which  strategy  to  chose  depends  on  the  data  at  hand  and  the  respective 
application.  If  the  data  quality  is  good  and  error  estimation  is  an  important  aspect,  the  over-all  strategy  is 
recommended.  If  the  data  quality  is  very  poor,  the  accumulation  strategy  can  be  the  only  way  to  obtain  a 
solution  at  all. 

3.2  Receiver  Autonomous  Integrity  Monitoring  (RAIM) 

If  the  signal  tracking  is  subject  to  a  temporary  discontinuity,  the  integer  part  of  the  ambiguity  N'  may 
change  while  its  fractional  part  remains  constant  -  a  so  called  cycle  slip  occurs.  Such  slip  biases  all 
subsequent  observations  of  the  respective  channel.  When  using  the  over-all  solution  strategy,  cycle  slips 
are  to  be  distinguished  from  other  outliers  and  excluded  for  the  present  epoch  and  the  future.  When  using 
the  accumulation  strategy,  both  outliers  and  cycle  slips  are  to  be  excluded  for  the  current  epoch  only  but 
discrimination  is  obsolete.  The  detection  and  exclusion  task  is  addressed  by  classical  RAIM  methods  as 
initially  suggested  by  [3]  for  code-based  single-point  processing  in  aviation  applications.  This  approach  is 
based  on  a  consistency  check  of  the  over-determined  navigation  solution  and  uses  the  following  residual 
based  test  statistic: 


ym  fjl 

RMS(f) «  ^ -  (35) 

\|  m  —  l 

fj  =n-\,n  DQj  _  «-!,»£,(£  / 

For  sampling  rates  of  1  Hz  or  higher,  RMS(f)  is  very  low  as  all  errors  other  than  the  ones  caused  by  cycle 
slips  or  other  outliers  virtually  vanish.  This  causes  distinct  spikes  in  the  test  statistic  in  the  presence  of 
outliers  allowing  for  their  detection  and  exclusion.  Note  that  this  approach  requires  the  incremental 
solution  which  therefore  has  to  be  calculated  even  when  applying  the  over-all  solution  strategy.  An 
elaborate  explanation  of  the  RAIM  method  in  time-differential  processing  is  beyond  the  scope  of  this 
tutorial.  The  interested  reader  is  referred  to  [11]. 

3.3  External  correction  data 

Also  non-outlying,  healthy  measurements  are  afflicted  by  errors  Djtotal .  When  working  with  the  present 
approach,  these  errors  are  growing  with  increasing  processing  time  spans.  For  attenuating  the  error  drift 
different  measures  are  taken  when  implementing  the  navigation  equations:  Tropospheric  signal  delays  are 
compensated  for  by  the  UNB3  model  given  in  [4],  Ionospheric  propagation  delays  (advances  when 
working  with  phase  observations)  are  accounted  for  by  a  thin  layer  model  as  indicated  in  the  standard 
literature,  e.g.  [9].  Here,  the  total  electron  count  (TEC)  is  computed  from  ionospheric  correction  maps 
provided  free  of  charge  by  the  Crustal  Dynamics  Data  Information  System  in  the  IONEX  format 
(www.cddis.gsfc.nasa.gov).  Ephemeris  error  reduction  is  realized  using  final  ephemeris  products 
published  online  by  the  International  GNSS  Service  (IGS,  http://igscb.jpl.nasa.gov).  Attenuating  the  range 
error  resulting  from  the  satellite  clock  frequency  offset,  which  integrates  to  the  satellite  clock  error,  is 
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achieved  by  the  use  of  30  s  sampled  clock  correction  data  published  by  IGS.  As  the  clock  error  is  growing 
faster  than  the  ephemeris  error,  these  corrections  are  particularly  important  when  working  with  time 
differences.  More  details  about  the  effect  of  the  described  correction  models  on  the  time  difference 
solution  is  given  by  [15].  With  Eq.  (29)  one  can  see  that  the  error  bl  V/geo  stemming  from  a  bias  in  the 

base  solution  becomes  maximum  if  A^is  collinear  with  e, .  According  to  [17]  |e;||  does  not  exceed 

0.00019  Is.  For  a  processing  interval  of  200  s  and  a  base  position  bias  of  4  m  the  resulting  range  error 
would  be  15  cm  in  this  worst  case  scenario.  Hence  it  is  important  to  apply  the  range  corrections  used  for 
time  differential  processing  to  attenuate  error  drift  also  when  calculating  the  base  position  via  code  based 
single  point  processing  for  absolute  error  limitation. 


4.0  PRACTICAL  VALIDATION 

4.1  Static  Experiment 

The  results  of  a  static  experiment  are  discussed  in  order  to  demonstrate  the  precision  achievable  by  the 
time  difference  approach  under  good  conditions.  Moreover,  the  effect  of  the  geometric  error  stemming 
from  an  offset  in  the  initial  base  position  is  analyzed.  The  test  was  performed  on  June  19th,  2006  on  a 
sports  field  in  the  north  of  Munich  with  scattered  buildings  nearby  (48°16’52.3”  N;  1 1°40’30.2”  E).  The 
low  cost  receiver  used  is  a  u-Blox  TIM-LP  single  frequency  module  integrated  in  an  evaluation  kit 
provided  by  the  manufacturer.  Zero  baseline  tests  show  the  standard  deviation  of  the  phase  measurements 
of  the  receiver  to  be  as  low  as  0.87  mm  [12],  For  processing,  ionospheric  correction  maps  (igsgO  1 90.06i), 
precise  ephemeris  (COD13584.EPH)  and  clock  corrections  (COD13584.CLK)  are  used.  A  time  interval  of 
5  min  is  processed,  which  is  significantly  longer  than  the  typical  intervals  to  be  analyzed  for  flight 
measurements. 

The  bold  red  line  of  the  upper  plot  of  Figure  6  shows  the  error  drift  of  the  over-all  solution  during  this 
time  interval,  i.e.  the  offset  of  the  static  relative  solution  from  zero.  Good  PDOP  vales  of  about  2.1,  see 
lower  plot,  limit  the  error  to  below  3  dm.  Besides  the  error,  the  estimate  of  the  error  is  illustrated  by  the 
upper  plot  of  Figure  6:  The  lowest,  dashed  green  line  is  the  root  mean  square  of  the  phase  range  residuals, 
the  light  blue  line  right  above  is  the  estimate  of  the  standard  deviation  of  the  measurements,  compare  Eq. 
(32).  Scaling  with  PDOP  yields  the  3D  position  error  estimate  as  given  by  Eq.  (34).  Considering  the 
simplifications  made  for  calculating  this  estimate  indicated  by  the  bold,  dotted,  grey  line,  it  matches  the 
actual  error  well.  However,  one  has  to  be  aware  that  the  quality  of  the  error  estimate  depends  on  the 
number  of  used  satellites  and  will  be  meaningless  for  the  case  of  only  four  used  satellites.  The  impact  of  a 
bias  in  the  base  position  was  discussed  theoretically  by  Eqs.  (25)  through  (29)  and  illustrated  by  Figure  4. 
This  effect  is  validated  by  adding  an  offset  vector  d\h  with  random  direction  but  defined  length  to  the  best 
estimate  of  the  initial  solution  as  used  above.  The  results  are  shown  by  Figure  7: 
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Figure  6:  Static  test:  3D  error  with  quality  assessment. 
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Figure  7.  Static  test:  impact  of  a  bias  in  the  base  position  on  relative  precision 


One  can  clearly  see  a  (linear)  relationship  between  the  offset  in  the  base  position  and  the  growing  drift  in 
the  time  difference  solution.  These  results  coincide  well  with  Eq.  (29)  and  underline  the  need  for  a  good 
estimate  of  the  base  position  and  time  solution.  However,  the  results  of  this  static  test  also  show  that  the 
precision  of  the  time  difference  solution  obtained  when  applying  all  mentioned  corrections  stays  in  the  low 
decimeter  (or  even  centimeter)  range  depending  on  the  time  spans  to  be  analyzed. 

4.2  Dynamic  Flight  Experiment 

HARDWARE.  The  test  was  performed  using  the  miniaturized  GPS  data-logger  depicted  in  Figure  8. 
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Figure  8:  Miniaturized  GPS  datalogger  with  25x25  mm2  patch  antenna  attached 

This  receiver  features  the  u-Blox  module  LEA-4T  capable  to  provide  raw  GPS  data.  Signals  are  acquired 
by  a  25  x  25  x4  mm3  passive  patch  antenna  placed  on  35  x  35  mm2  ground  plate.  In  addition  to  the 
onboard  PT  solution  available  with  up  to  4  Hz,  non  processed  raw  data  with  a  sampling  rate  of  up  to  10  Hz 
can  be  logged  within  the  8  MB  internal  flash  memory  (approximately  40  min  of  recording  at  10  Hz).  The 
quality  of  the  raw  data  is  expected  to  be  comparable  to  those  of  the  TIM-LP  receiver  used  for  the  static 
test.  The  form  factor  of  the  logger  is  as  small  as  44x21x4  mnT  at  a  weight  of  4.25  g  only  without  battery. 
Adding  the  patch  antenna  (14.15  g)  and  a  320  mAh  3.7  V  lithium-polymeric  battery  (7.1  g)  increased 
weight  of  the  construction  up  to  25.5  g.  The  power  consumption  is  approximately  40  mA.  The  available 
logging  memory  is  extended  by  connecting  an  external  serial  datalogger  featuring  a  2  GB  SD  memory 
card.  The  volume  of  the  serial  data  logger  is  50  x  36  x5  mm3,  its  weight  is  6  g  and  its  power  consumption 
is  6.7  mA.  2  GB  of  memory  is  sufficient  to  store  10  Hz  raw  GPS  data  during  8  days  approximately. 

FLIGHT  TEST  WITH  ACROBATIC  PLANE.  The  flight  was  performed  by  the  acrobatic  plane  Mii30 
“Schlacro”  developed  and  built  by  the  academic  aviation  association  AKAFLIEG  of  TU  Miinchen.  The 
receiver  was  mounted  on  the  outer  left  wing  tip  of  the  aircraft,  as  illustrated  by  Figure  9. 


Figure  9:  Mu30  flight  test:  Aircraft  Mu30  “Schlacro” 
with  GPS  logger  mounted  on  left  wing  tip 


The  10  Hz  phase  data  collected  during  a  S-sequence  maneuver  flown  close  to  the  airfield 
Fiirstenfeldbruck,  Gennany,  on  August  15th,  2008,  is  analyzed.  The  2D  trajectory  and  the  vertical  profile 
of  the  maneuver  are  given  by  Figure  10. 
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(a)  2D  trajectory.  (1)  Start  of  precise 

relative  processing;  (2)  end  of  precise  (b)  Vertical  profile  relative  to  the 

relative  processing  starting  point  (1)  of  the  maneuver. 

Total  time  history  (to  the  left)  and 
zoomed  section 


Figure  10:  Mu30  flight  test:  flight  trajectory  during  S-sequence  maneuver. 


When  processing,  precise  ephemeris,  high  rate  clock  solutions  and  ionospheric  correction  maps  are  used. 
The  quality  of  the  collected  phase  measurements  is  good,  the  root  mean  square  of  the  C/NO  values  yields 
48.1  dBHz.  As  a  consequence,  the  residual  level  and  the  estimate  of  the  measurement  standard  deviation 
dropping  out  from  the  least  squares  solver  are  very  low,  see  light  green  and  blue  lines  in  the  upper  plot  of 
Figure  11. 


tau  [s] 


Figure  11:  Mu30  flight  test:  quality  analysis  of  the  time  difference  solution 


The  geometry  of  the  available  satellite  subset  results  in  PDOP  values  of  up  to  4.4,  lower  plot  of  the  figure. 
These  values  are  combined  to  the  3D  error  estimate,  Eq.  (34),  indicated  by  the  bold  grey  line  of  the  upper 
plot  of  Figure  11  The  precision  is  estimated  to  be  better  than  17  cm.  For  validation  purposes,  a  second 
receiver  (u-Blox  TIM-LP)  was  mounted  close  to  the  airfield  serving  as  a  base  station  for  RTK  processing 
and  static  initialization  was  provided.  The  fixed  solution  is  computed  in  a  differential  (RTK)  mode  using 
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the  data  of  the  base  receiver.  The  double-differenced  phase  and  code  data  of  both  receivers  are  processed 
using  inhouse  Kalman  filter  software  of  the  Department  of  Earth  Observation  &  Space  Systems  of  Delft 
University  of  Technology  [14].  It  is  well  worth  to  note  that  aside  from  serving  as  a  validation  of  the  time 
difference  method,  the  RTK  solution  also  demonstrates  successful  ambiguity  resolution  when  processing 
data  of  truly  low  cost  receivers  collected  in  highly  dynamic  applications.  The  deviation  of  the  obtained 
results  from  the  time  difference  solution  is  shown  by  Figure  12.  The  relative  error  components  in  local 
tangent  coordinates  do  not  exhibit  a  special  trend  to  the  height  direction  and  the  absolute  value  of  the  3D 
relative  offset  stays  below  13  cm.  This  result  confirms  the  high  precision  achievable  by  time  differences. 


Figure  12:  Mu30  flight  test:  comparison  of  the  time  difference  method 
with  a  RTK  solution  generated  by  DEOS  of  TU  Delft  [14] 


5.0  IN  FLIGHT  MEASUREMENTS  OF  DYNAMIC  SOARING  OF 
ALBATROSSES 

The  Wandering  Albatross  lives  predominantly  at  very  remote  places  of  the  southern  hemisphere.  As  a 
consequence,  the  realization  of  our  measurement  campaign  during  the  arctic  summer  2008/09  at 
Kerguelen  Archipelago  (Figure  13  (a))  was  only  possible  in  cooperation  with  our  partners  from  the  Centre 
d’Ecologie  Fonctionnelle  et  Evolutive  of  the  French  Centre  National  de  la  Recherche  Scientifique 
contributing  the  indispensable  ornithological  and  logistical  expertise.  As  recovering  the  receivers  after 
application  on  the  birds  was  required,  measuring  took  place  during  the  breeding  period,  when  the  birds 
take  turns  leaving  their  nest  for  up  to  30  day  lasting  foraging  trips  over  open  waters. 
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(a)  Kerguelen  (b)  Wandering  Albatross  equipped  with  GPS  receiver 

Archipelago 


Figure  13:  GPS  measurement  campaign  during  the  arctic  summer  2008/09 

The  birds  where  equipped  with  hardware  comparable  to  the  receiver  presented  in  section  0.  Power  was 
supplied  by  three  lithium-thionyl-Chloride  AA-type  primary  cells.  The  mass  of  the  complete  logging  unit 
did  not  exceed  100  g.  Some  receivers  additionally  featured  3  axis  MEMS  accelerometers.  Mounting  on  the 
bird  was  realized  using  adhesive  tape,  see  Figure  13(b).  16  flights  /  flight  sections  where  successfully 
recorded,  two  of  them  are  shown  in  Figure  14  (a)  and  (b). 
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(a)  3  day  foraging  trip 


(b)  First  6  days  of  a  30  day  trip 


(c)  3D  view  of  the  bird’s  trajectory  with  dynamic  soaring  cycles 


Figure  14:  GPS  tracks  of  foraging  Wandering  Albatrosses 


5.1  Analysis  of  an  Exemplary  Dynamic  Soaring  Cycle 

For  understanding  the  dynamic  soaring  technique,  the  global  trajectories  are  of  secondary  interest  only  but 
individual  flight  cycles  as  shown  in  Figure  14  (c)  are  to  be  analyzed  precisely  using  the  time  differential 
technique.  The  processing  results  of  an  exemplary  cycle  are  discussed  in  the  below. 

Figure  15  shows  a  randomly  chosen  flight  cycle.  In  the  left-hand  side  plot,  the  vertical  component  is 
scaled  by  a  factor  of  three  for  improving  the  perspective  perception.  The  time  history  of  the  vertical 
component  of  the  16  second  lasting  maneuver  already  gives  an  impression  of  the  precision  of  the  time- 
differential  solution,  Figure  15  (b).  The  solution  was  calculated  applying  the  accumulation  strategy;  5  to  9 
satellites  where  available  yielding  PDOP  values  ranging  from  1.8  to  5.5.  The  error  drift  is  estimated  to 
reach  up  to  0.9  dm  -  an  estimate  which  is  probably  overly  pessimistic,  as  mentioned  below. 
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(a)  3D  trajectory  of  an  individual 
dynamic  soaring  cycle 


(b)  Time  history  of  the  vertical 
component 


Figure  15:  Dynamic  soaring  cycle  processed  using  the  time-differential  technique 


For  further  analysis  the  position  fixes  are  transformed  to  an  analytical  form  using  cubic  smoothing  splines 
as  a  first  approach.  Such  interpolation  is  shown  for  the  vertical  component  by  Figure  16.  Flere  the 
smoothing  parameter  is  chosen  in  a  way  that  the  spline  residuals  (bottom  plot)  are  of  the  same  order  of 
magnitude  than  the  3D  estimate  of  the  noise  component  calculated  according  to  Eq.  (34)  using  inter-epoch 
residuals,  center  plot.  Staying  below  2  cm,  the  residual  level  is  very  low  considering  the  high  dynamics  of 
the  maneuver. 


1 

•  original  data  / 

- smoothing  spline  fit  / 

cl _ I _ I _ I _ 

0  5  10  15 


tau  [m] 


Figure  16:  Interpolation  and  smoothing 


The  fitted  position  solution  can  now  be  used  for  velocity  calculation  -  the  blue  line  of  the  upper  plot  of 
Figure  17  illustrates  the  speed  relative  to  the  ground  resulting  from  analytical  derivation.  The  red  dots  of 
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the  same  plot  are  obtained  by  directly  differentiating  the  position  data  (using  MATLAB’s  diff  command) 
without  precedent  smoothing.  (Note  that  both  solutions  are  significantly  smoother  than  the  solution 
obtained  using  the  recorded  Doppler  measurements  for  single  point  processing,  illustrated  by  the  black 
circles).  The  lower  plot  shows  the  corresponding  absolute  acceleration:  again  the  red  dots  are  the  result  of 
MATLAB’s  diff  command  applied  twice  to  the  position  data  whereas  the  blue  line  is  the  second  derivative 
of  the  analytical  fitting  results.  The  green  line  is  the  absolute  acceleration  as  measured  by  the  3D  MEMS 
accelerometer,  also  available  with  10  Hz.  Contrarily  to  the  acceleration  derived  from  the  GPS  position, 
this  value  comprises  the  gravity  vector.  The  green  and  the  blue  line  are  not  parallel  due  to  the  tilting  of  the 
accelerometer  with  respect  to  the  east-north-up  coordinate  system.  For  the  time  being,  the  accelerometer 
results  are  only  shown  here  as  a  coarse  consistency  check. 


Figure  17:  Dynamic  soaring  cycle  processed  using  the  time-differential  technique 


5.2  Outlook 

This  exemplary  analysis  aims  to  demonstrate  the  potential  of  the  time-difference  method  applied  to  data 
gathered  by  miniaturized  GPS  receivers  under  rather  difficult  field  conditions.  Future  work  comprises 
aspects  such  as  the  development  of  more  sophisticated  interpolation  methods,  the  recalibration  and 
analysis  of  the  used  MEMS  accelerometers  (little  is  known  about  the  quality  of  theses  sensors  by  now) 
and  of  course  the  pursuit  of  the  flight  mechanical  evaluation  of  the  measured  data. 


6.0  CONCLUSIONS 

Carrier  phase  time  differences  are  used  in  a  diversity  of  applications  and  recently  also  gain  importance  in 
safety  critical  applications  such  as  receiver  autonomous  integrity  monitoring  for  aviation  applications.  A 
method  to  utilize  the  fact  of  canceling  ambiguities  when  forming  time-differences  as  a  standalone 
possibility  to  precisely  process  kinematic  measurements  over  time  spans  of  several  minutes  has  been 
presented  within  this  paper. 
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The  navigation  equations  underlying  this  approach  were  derived  and  theoretical  aspects  of  various  error 
sources  were  reflected  including  an  estimation  of  their  impact  on  the  final  base  vector  precision.  Further  a 
short  introduction  to  integrity  monitoring  and  the  outline  of  two  different  solution  strategies  for  trajectory 
measurement  was  given. 

The  practical  validation  of  static  data  gathered  under  favorable  conditions  confirmed  that  the  temporal 
correlation  of  the  remaining  range  errors  is  high.  This  allowed  for  processing  time  spans  of  up  to  5  min 
with  relative  aberrations  below  30  cm.  Moreover,  the  test  showed  that  the  derived  error  estimate 
principally  is  a  suitable  means  to  assess  the  quality  of  the  navigation  solution.  The  impact  of  position 
biases  of  the  base  epoch  was  also  demonstrated  with  the  static  data.  Such  offsets  cause  an  additional 
“geometric”  range  error  which  can  significantly  distort  the  baseline.  Applying  the  correction  data  used  to 
attenuate  the  error  drift  also  when  calculating  the  initial  position  via  code  based  single  point  processing  is 
a  suitable  means  to  overcome  this  problem.  Dynamic  flight  test  results  performed  with  a  completely  self- 
contained,  miniaturized  GPS  logging  device  featuring  a  massmarket  LI  single  frequency  GPS  module 
were  presented.  Comparison  with  a  fixed  ambiguity  RTK  solution  demonstrated  that  decimeter  precision 
can  be  achieved  even  under  very  dynamic  conditions. 

The  development  of  the  time-difference  method  and  the  according  hardware  was  triggered  by  the 
objective  to  measure  the  dynamic  soaring  of  Albatrosses  using  GPS.  This  project  was  successfully  realized 
during  winter  2008/09.  First  results  from  this  measurement  campaign  have  been  presented  in  this  paper. 
The  raw-data  evaluation  based  on  the  time-difference  method  was  demonstrated  under  navigational 
aspects  and  the  subsequent  data  treatment  yielding  speed  and  even  acceleration  was  shown. 

The  time-difference  method  is  both  unconventional  and  easy  to  apply  in  the  field:  For  a  user  who  is  well 
aware  of  the  pitfalls  of  raw  phase  data  processing,  neither  a  second,  nearby  base  station  along  with  (static) 
initialization  procedures  nor  high  quality  receiver  technology  is  required  to  achieve  a  truly  phase  based, 
stand-alone  navigation  solution  precise  to  the  low  decimeter  range.  These  advantages  propose  the 
application  of  the  present  approach  not  only  for  the  unconventional  task  of  the  flight  mechanical 
investigation  of  birds  but  for  a  wide  panoply  of  other  technical  applications.  One  might  think  of  measuring 
take-off  or  landing  distances  of  aircraft  during  the  certifications  process  or  the  quality  assessment  of 
trajectories  autonomously  flown  by  UAVs. 
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